{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "4690e96d-4c34-42f1-aad0-3d54e3939d10",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of nonzero wavedecn coefficients = 535\n",
      "Number of nonzero fswavedecn coefficients = 280\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_48705/3245871848.py:48: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.\n",
      "  img[slices] = val\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWIAAACeCAYAAAAIewCsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAVMklEQVR4nO3debhcRZ3G8e/LjiQQ9iVi2JFN4kKCQAIqyrBnBpVghKCC6DxG9GFRGISoQEQUGUZ5cIQBZWdgWJTgREVIQDYFUUEUBGIwLCEJkLAM22/+qOrk5OQu3XdJ9b28n+c5z72nqs451V3dv66z1VFEYGZm5SxXugJmZm91DsRmZoU5EJuZFeZAbGZWmAOxmVlhDsRmZoU5EJuZFeZAbC2R9LikZyStVkk7QtItBavVK5IukvSqpAV5+pOkKZLWyPkrSFooaXRlmQmSooO0hySdl8svzOt9rTJ/U4nX2JlB3p4LK9PBOe8zuY0WSHpa0lRJQyW9P6ctX1nPjzpJO0/SA5V1vyHplcr8ia3W2YHYemJ54OjSleiMpBV6sNi3I2IosC7wKWBn4HZJq0XE68AdwNhK+bHAQx2kTY+Iz0XEkIgYApwOXNmYj4i9e/Ka+tlgbc8hlelKSbuT2uOQ3NbbAFfm8r8lxcP3VNYxBniiltZo4+0qbTwD+EJlW6e3WlkHYuuJM4FjJQ3rKFPSLpLukfR8/rtLJe8WSd+UdHvubUyTtE7O+36tF/O6pMk5byNJ10iaI+kxSV+srHOypKslXSLpBeDwXP4GSfMkPSLpyGZeWES8EhH3AAcAa5OCMsB0lgy6Y4AzOkib3sx22sygbc+anYA7IuI+gIiYFxE/jogFEfEacCe5PSWtB6wEXFVL24p+aGMHYuuJ3wK3AMfWMyStBdwInEMKZGcBN0pau1LsE6QA1/iwHwsQEV+o9DJ2A+YD10taDvgpcD8wHPgQ8CVJe1XWeSBwNTAMuBS4gtSb2Qj4KHC6pA82+wIjYgHwC1JwhfTl21XScjnQrEb6ko6qpG3DwAzEg749s7uAvSR9XdKuklau5Vd/bMcCt+WpmvZYRDzR4na75UBsPXUyMEnSurX0fYGHI+LiiHg9Ii4n7cLvXylzYUT8NSJeJgWzkdUV5HVeB0zKvZedgHUj4hsR8WpEPAr8CBhfWeyOiLguIt4E1gF2Bb6Se7i/B84HDmvxNc4G1sr/3wW8DdiBFJxvi4iXgMcqaY9HxN9b3Ea7GGzteayk5/L0LEBEzAD+hXSo4UZgrqSzKseAbwV2kyRSe84gHZLauZJ2axfb7DEHYuuRiPgT8DPgq7WsjYCZtbSZpJ5Pw1OV/18ChjRmJK1I6gldFhFX5OQRwEaVL9ZzwInA+pX1zKrVYV7u1XZWh2YMB+ZBOmQB3E3qFY0lfUlhcY9pLAOzNwwMyvb8TkQMy9M6ldd5U0TsT/qBPRA4HDgiZ9+Z6749uY0jYmGuSyOtX9rYgdh64xTgSJb8QswmfdGq3gH8o8l1/gfwAnBSJW0WaZdwWGUaGhH7VMpUhxGcDawlaWgP64CkIcCeLA64sHjXdUwlfUYlbcAG4mzQtmddRLwZEb8CbiYF2caP7T2k3v6GEfFQLj4jp70LB2JrNxHxCOms8xcryVOBrSR9Qumyr4OBbUm9rS5JOgrYHZiQd0kb7gYWSPqKpFUlLS9pe0k7dVKvWcBvgCmSVpH0LuAzwCVN1GFlSe8l7UrPBy6sZE8HPgBsDDyY024H9iDtjg/oQDwY27NWnwMljZe0ppJRuX53VopNJ11B8ptK2m057cmI+Fsr22yWA7H11jdIJ64AiIi5wH7AMcBc4Hhgv4h4tol1HQJsBsyunGk/MSLeyOscSTom+yzpGOEa3axrE1Jv6lrglIj4ZRflj5e0INf5J8DvgF0i4sVKmd/kbd4VeSDv/LrmAM9ExMNNvMZ2N1jasyPzST3+h0m99EuAMyPi0kqZW0knHW+rpN2W06p7R31KHhjezKws94jNzApzIDYzK8yB2MysMAdiM7PCHIitbUlaX9L0PIbBd/MlRxdKmi/pbkljJP2lifVMkDRtWdTZuud2XZqvmrC2JelrwLuBgyIiJI0BLge2rl1WtizrFMCW+Zpb6wG369LcI7Z2NgJ4MBb3FkaQxnMo8mW1PuN2rYsIT576ZCLdcfY/pBsc5gLfJ/3Yn0QaG+AZ0s0Sa1SW2Zl0o8RzpNG49sjpFwGvAa8CC4GjgFeAN/L810l3tD3R1fZz+uGkQXoa5d5JGlltHvAX4OOVvIuAH5AGhVlAGuxn85w3nXTr7Yu5DgeTBqT5Wa7/PNJF/8uVbgu368Bq1+KN7GlwTKTBxe8Hvke6M2sV0tCHnwYeId1hNSR/oS7OywzPX6x98hf7w3l+3Zx/EXBqZRv1L96iL2xn268vl/NmkYZtXIG0i/wssG1lm3OBUTn/UuCKyjYD2KIyPwU4D1gxT2PIh/wGw+R2XTbt6kMT1ldGkUbJOi4iXow0XOFtwATgrIh4NNJIVicA45WeuvBJYGpETI00CMsvSGPj7tPZRnqw/br9SLvBF0Ya1vE+4BrgY5Uy10bE3ZGezHEptWEda14DNgRGRMRrETEj8jd5kHC7LoN2dSC2vrIxMDN/yKvqwyjOJPVI1icdG/xYbTjE3UhfgL7aft0IYHRtmxOADSplOh3WsQNnknqG0yQ9Kqk+jORA53ZdBu3ak2dBmXVkFvAOSSvUvjT1YRTfAbwOPJ2XuTgievLYm2a331G5WyPiw32wTSKNkXsMcIyk7YGbJd0TaYjFwcDtugza1T1i6yt3A08C35K0Wh6ucFfSZUlflrRpHuO38TDN10mjX+0vaa88FOIqkvaQ9PY+3H7dz0jDOh4qacU87SRpmya38zTpuCgAkvaTtIUkAc+TTjq92dnCA5DbdRm0qwOx9YlIQxvuD2wB/J30fLGDgf8CLiadmX6MdIZ8Ul5mFukpCSeSzojPAo6jB5/LLrZfL7cA+AjpsTyzSburZwD155d1ZjLw47z7+3FgS+CXpLPtdwDnRsSvW61/u3K7Lpt29Q0dZmaFuUdsZlaYA7GZWWEOxGZmhTkQm5kV1u+BOD8wcLPuS5Yj6bw8IpT1M0lTJH2pdD1aIekaSXuXrkc7G6DtOknSGaXrATQ/1gRwC+kpqCv31/3WXWx7ifvA22UCHgf2LF2PgTIB6wL/AFbN83vktj23Vu424PAW170H6TrPhZVpYiX/EtL1qC8AfwWOqORtkutRXfZrlfxRwO9Kv3/tOvVzu+6bl3uOdEna+cDQSv5awJWkcSSeJd26vHqtbX9NupPuoer3lTRuxRPAeqXfw6Z6xJI2IQ16EcABzSyzLOX72639HU4ag+DlStqLwKH5M9ZbsyNiSGX6cSVvCrBJRKxO+gyfKum9teWHVZb9ZiMxIu4GVpf0vj6o42B0OP3XrmsAp5Juqd6GNKDQmZX8U4E1gU2BzUm3WE+u5F8O3AesDfwbcLWkdQEi4hXgJuCwXtax15o9NHEYcCdpBKOJ1QxJF0n6gaQb84j7d0navJIfkrbI/+8j6cFc7h+Sjq2UO1LSI5LmSbpB0kY5fXoucn8+zHFwvkvnCUlfkfQUcKGklSWdLWl2ns6WtHJeR6P8iZKelfS4pAm113BqZf5ASb+X9IKkv0n6p+7eIEmHS7pd0vfyReGPStolp8+S9IykiZXy+0q6L29jlqTJtfUdJmmmpLmSvpbrvGfOW07SV3Pd5kq6StJa3dWxDewN3FpLe470uTqlPzccEQ9ExP81ZvO0eReL1N1C6p3Z0vqtXSPisoj4eUS8FBHzgR8B1TvrNgWui4gXIuJ54FpgOwBJWwHvAU6JiJcj4hrgj8BBleVvoQ3atZVAfGme9pK0fi1/PGkc0TVJA2Wc1sl6LgCOioihwPbAzQCSPkjqsXycNDDITOAKgIgYm5fdMfdUrszzG5B2S0YAnyX92u1MGlFpR9Lu5EmVbW9AGmN0OOnH5D8lbV2voKRRpLFVjwOGAWNJhyCaMRr4A+nX97L8GnYi3RX0SeD7+XZQSD2Gw/I29gU+L2lcrsO2wLmkQUs2JPUKhle2MwkYB+xO6inMJ4212u52II0TW3cacFAn7bGbKgO5dDDtVim+nqSnJT2WfxBXq63rXEmNXdQngam1zc3MP9gXSlqnlvdn0ufKltbf7Vo1FnigMv8DYD9Ja0pakxRkb8p52wGPRrrrruH+nN7QHu3axDGa3UhDwq2T5x8CvlzJvwg4vzK/D/BQZX7R8V3SLYpHUTmGk9MvAL5dmR+St7lJfR2VY1CvAqtU0v4G7FOZ34s0LF6j/OvAapX8q8jHAamMjwr8EPhek8evHicfcyLtnj1cydsh13v9StpcYGQn6zq7sV3gZODySt7b8uttbOvPwIcq+Rvm92uF0se6unm/XgPeWWvHxriz3yaNVQA9O5a4AbAtqXOxKenW2x92UG75/Jk+CVix8nl7H4tHD7sa+N/ackcCN5d+D9tx6s92rW3nw6ROx1aVtI1ItyK/madfACvlvEOBO2vrOA24qDK/JfBG6fewmR7xRGBaRDyb5y+jdniC5oeXO4gUqGdKulXS+3P6EkPqRRrfdC5L9gLr5kQ6xtPQ0bB8G1Xm58eSj2Kp5zdsTArqPfF05f+XASKinjYEQNJoSb+WNEfS88DnSD12cr1mNRaKiJdI70fDCOBaLR7u78+kQUnqeyrtZj4wtJO8M0h7Wz3qnUTEUxHxYKTxbx8DjmfJXdBGuTcijWf7duDzOW1hRPw20ji2TwNfAD4iqVrXoaTdbVtav7Vrg6SdSbHnoxHx10rWVaSTr0OB1Unf3Uty3sKcVrU66QkdDUNJg/oU1WUglrQq6XDB7pKeysdjvwzs2JM3NiLuiYgDgfWA60hvItSG1Mu7lGuTzsR2urrafEfD8s2uzK9Z21Wt5zfMorVjhz11GXADsHFErEF6GoBy3pOkQAEsaoe1a3XcOyKGVaZVIqKr96sd/AHYqqOMiJhL2iv4ZjVd6Ym+C7uYxnSyraDrz/cKdN7Ojc9WdfltSLu1trR+bVdJ7yZ9Vz4dSw9DOZK05/Ni7sCdx+IB6B8ANqv9oO7Ikoc22qJdu+sRjyP1tLYlveCRpIrPoMUzjZJWUnr89RoR8RrpMqLGsHKXA5+SNFLpBNvpwF0R8XjOX2KIuk5cDpwkad18fO9kFv8yNnw912MMaUT//+5gPRfkunwonxQbLumdrbzWJg0F5kXEK/m49CcqeVeThhHcRdJKpLPAquSfB5wmaQRAfs0H9kMd+9pU0nHtzpwF7EL6jAEQ6ckIQ7qYZgBI+oCkEUo2Br4FXJ/z1pM0XtIQpWEZ9wIOAX6V80dL2jq399rAOcAtkU7+NOzO4mOPtqT+bNftgZ8DkyLipx2s+x7gCEmr5g7LZ0k/DOSe8++BU5SGz/xn4F2kJ3c0tEW7dheIJwIXRsTf867fUxHxFOnhgRPU+mVjhwKPS3qBtCs+ASAifgl8jfQGPUnqqYyvLDeZJYeo68ippMex/IF0ZvTenNbwFGkXajbppOPnIuKh+koiXar0KdIzsp4nnQ0eUS/XB/4V+IakBaQfjcbeARHxAOmE3BWk92Mh6QGNjbP+/07qIUzLy99JOlHY7n4C7JO/MEuJiBdIxxR7cgXIu0kPq3wx//0j8MXGqkmHIZ4gfQa+A3wpIm7I+ZuRvuwLgD+R3udDGiuWtBOwMH82bGn92a7HkK5TvqDSW672aD9Nulb4CdIe9GYseeh0POn4/3zSj/NHI2IOgKRVSL3n6mWORbwlhsGUtAdwSUT0ZGDq4pSutHgO2DIf/xywJJ0OPBMRZ5euS7MkXQNcEBH1qywsG6DtOol0aPD44nVxIG5PkvYn7ToL+C6px/ueeCs0mNlbjAf9aV8Hkg6jzCZdYjPeQdhscHpL9IjNzNqZe8RmZoW1PFjOEUcc0XQXes8991wqbbvttuugZN8ZN25cS+V33TXdtr7VVh1eBrmE0aNHc/LJJ/ekWv3ijjvuUPelWtJh20oiIpDS5hp7UdX5jv7vcSWaWEdvttPKsn3xenoiIvpyo27XPtxOb3TWru4Rm5kV5kBsZlaYA7GZWWEOxGZmhTkQm5kV5kcM2SJdnUWu5pU422w919N27es2b2YdvdlOK8u222fYPWJbpLNBqxt59XI2MPS0XTv7vzf16IsyfbFsu31+3SO2RdwjHpzcI+7b7fQH94jNzApzIDYzK8yB2BbxMeLByceI+3Y7/cHHiG2RdjmWaH2rXdrVx4g75x6xmVlhDsRmZoU5EJuZFeZAbGZWmAOxmVlhDsRmZoU5EJuZFdbydcQ77rhjrza49dZb92r57owePbql8gcccAAAJ5xwQrdl33zzTebMmdOjepmZdcY9YjOzwhyIzcwKcyA2MyvMgdjMrLCWT9ZNmjSp6bJHH330Umljx45tdZMtGTZsWEvlx40bB8C0adO6LXvvvfdy3HHH9aBW1or+HhzGrN24R2xmVpgDsZlZYQ7EZmaFeWB4W2RZPLWgmWO7EdFtuWbKWNJZu0pa4n1slKvOd/R/b+rRn+3ayrLt9vlxILZF2uWD6ZN1fctP6Ojb7fQHH5owMyvMgdjMrDAHYjOzwhyIzcwKcyA2MyvMgdjMrDAHYjOzwhyIzcwKcyA2MyvMgdjMrDAHYjOzwhyIzcwKa3nQn1YeJz958uSl0q6//vpWN9mS888/v6Xyo0aNAuD+++/vtuw555zDyJEje1ItM7NOuUdsZlaYA7GZWWEOxGZmhTkQm5kV5kBsZlaYA7G95QzkR+rY4ORAbGZWmAOxmVlhDsRmZoU5EJuZFeZAbGZWmAOxmVlhDsRmZoU5EJuZFeZAbGZWWMvjEdvgFRGlq9BrzdwJFxFN3zHXStl21Vm7Slri9TXKVec7+r839ehuHb3ZzkBuVwdiW6SdPpj96a12i3NXr6GaVy/XVV5f16MvtjOQ29WHJszMCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8LU6lMZpkyZ0vQCw4cPXypt4sSJLW2vVWeddVZL5U877TQAzj777G7LTp06lRNOOKEn1eoXO+ywQ5+Obi1p4D+iY5CIiD5rW7dr++isXd0jNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrzIHYzKwwB2Izs8IciM3MCnMgNjMrrOWnOJuZWd9yj9jMrDAHYjOzwhyIzcwKcyA2MyvMgdjMrDAHYjOzwv4fvMaeKaDJoQoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"Using the FSWT to process anistropic images.\n",
    "\n",
    "In this demo, an anisotropic piecewise-constant image is transformed by the\n",
    "standard DWT and the fully-separable DWT. The 'Haar' wavelet gives a sparse\n",
    "representation for such piecewise constant signals (detail coefficients are\n",
    "only non-zero near edges).\n",
    "\n",
    "For such anistropic signals, the number of non-zero coefficients will be lower\n",
    "for the fully separable DWT than for the isotropic one.\n",
    "\n",
    "This example is inspired by the following publication where it is proven that\n",
    "the FSWT gives a sparser representation than the DWT for this class of\n",
    "anistropic images:\n",
    "\n",
    ".. V Velisavljevic, B Beferull-Lozano, M Vetterli and PL Dragotti.\n",
    "   Directionlets: Anisotropic Multidirectional Representation With\n",
    "   Separable Filtering. IEEE Transactions on Image Processing, Vol. 15,\n",
    "   No. 7, July 2006.\n",
    "\n",
    "\"\"\"\n",
    "\n",
    "import numpy as np\n",
    "import pywt\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "\n",
    "def mondrian(shape=(256, 256), nx=5, ny=8, seed=4):\n",
    "    \"\"\" Piecewise-constant image (reminiscent of Dutch painter Piet Mondrian's\n",
    "    geometrical period).\n",
    "    \"\"\"\n",
    "    rstate = np.random.RandomState(seed)\n",
    "    min_dx = 0\n",
    "    while(min_dx < 3):\n",
    "        xp = np.sort(np.round(rstate.rand(nx-1)*shape[0]).astype(np.int64))\n",
    "        xp = np.concatenate(((0, ), xp, (shape[0], )))\n",
    "        min_dx = np.min(np.diff(xp))\n",
    "    min_dy = 0\n",
    "    while(min_dy < 3):\n",
    "        yp = np.sort(np.round(rstate.rand(ny-1)*shape[1]).astype(np.int64))\n",
    "        yp = np.concatenate(((0, ), yp, (shape[1], )))\n",
    "        min_dy = np.min(np.diff(yp))\n",
    "    img = np.zeros(shape)\n",
    "    for ix, x in enumerate(xp[:-1]):\n",
    "        for iy, y in enumerate(yp[:-1]):\n",
    "            slices = [slice(x, xp[ix+1]), slice(y, yp[iy+1])]\n",
    "            val = rstate.rand(1)[0]\n",
    "            img[slices] = val\n",
    "    return img\n",
    "\n",
    "\n",
    "# create an anisotropic piecewise constant image\n",
    "img = mondrian((128, 128))\n",
    "\n",
    "# perform DWT\n",
    "coeffs_dwt = pywt.wavedecn(img, wavelet='db1', level=None)\n",
    "\n",
    "# convert coefficient dictionary to a single array\n",
    "coeff_array_dwt, _ = pywt.coeffs_to_array(coeffs_dwt)\n",
    "\n",
    "# perform fully seperable DWT\n",
    "fswavedecn_result = pywt.fswavedecn(img, wavelet='db1')\n",
    "\n",
    "nnz_dwt = np.sum(coeff_array_dwt != 0)\n",
    "nnz_fswavedecn = np.sum(fswavedecn_result.coeffs != 0)\n",
    "\n",
    "print(\"Number of nonzero wavedecn coefficients = {}\".format(np.sum(nnz_dwt)))\n",
    "print(\"Number of nonzero fswavedecn coefficients = {}\".format(np.sum(nnz_fswavedecn)))\n",
    "\n",
    "img = mondrian()\n",
    "fig, axes = plt.subplots(1, 3)\n",
    "imshow_kwargs = dict(cmap=plt.cm.gray, interpolation='nearest')\n",
    "axes[0].imshow(img, **imshow_kwargs)\n",
    "axes[0].set_title('Anisotropic Image')\n",
    "axes[1].imshow(coeff_array_dwt != 0, **imshow_kwargs)\n",
    "axes[1].set_title('Nonzero DWT\\ncoefficients\\n(N={})'.format(nnz_dwt))\n",
    "axes[2].imshow(fswavedecn_result.coeffs != 0, **imshow_kwargs)\n",
    "axes[2].set_title('Nonzero FSWT\\ncoefficients\\n(N={})'.format(nnz_fswavedecn))\n",
    "for ax in axes:\n",
    "    ax.set_axis_off()\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb6e8b72-fe4b-4c94-a3b7-5108ecec34a9",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
